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Abstract 



Wc developed a new direct-tree hybrid iV-body algorithm for fully self-consistent iV-body simulations 
' (— I '' of star clusters in their parent galaxies. In such simulations, star clusters need high accuracy, while galaxies 

^ |. need a fast scheme because of the large number of the particles required to model it. In our new algorithm, 

the internal motion of the star cluster is calculated accurately using the direct Hermite scheme with 
individual timesteps and all other motions are calculated using the tree code with second-order leapfrog 
^ ' integrator. The direct and tree schemes are combined using an extension of the mixed variable symplectic 

Ci . (MVS) scheme. Thus, the Hamiltonian corresponding to everything other than the internal motion of 

the star cluster is integrated with the leapfrog, which is symplectic. Using this algorithm, we performed 
fully self-consistent A^-body simulations of star clusters in their parent galaxy. The internal and orbital 
^ ■ evolutions of the star cluster agreed well with those obtained using the direct scheme. We also performed 

' fully self-consistent A^-body simulation for large- iV models {N = 2 x 10®). In this case, the calculation 

lO ■ speed was seven times faster than what would be if the direct scheme was used. 

' Key words: stellar dynamics — methods: numerical, n-body simulations — galaxies: star clusters 

(N 

■ 1. Introduction respect to its own gravity and starts to fragment. The 
f — I fragments collapse and star formation occurs. However, 

. Very young and massive stars have been found in the this model is problematic. Observations have shown that 
central parsec of the Galaxy (Krabbe et al. 1995; Paumard two disks are at large angles with respect to each other and 
et al. 2001; Paumard et al. 2006). These stars are very these stars on the disks formed almost at the same time. 

■ young (< 10 Myr) (Paumard et al. 2001; Ghcz 2003) and In this in situ formation scenario, two disks at large angles 
JH [ lie on two disks (Paumard et al. 2006). The disks rotate must have existed within 1 Myr, unless one disk changed 

around the central black hole (BH) and are at a large an- its orientation through, for example, external perturba- 
gle with each other. One disk rotates clockwise in projec- tion. 

tion, the other counterclockwise, although the existance of The star cluster inspiral model was proposed by 
the counterclockwise disk is controversial (Lu et al. 2006). Gerhard (2001). A star cluster is formed at a distance 
These disks are coeval to within 1 Myr. of tens of parsecs from the Galactic center (GC) and spi- 

Howcver, in situ formation of these stars seems unlikely raled into the GC due to dynamical friction before be- 
since the strong tidal field of the central super massive ing disrupted by the tidal field of the central black hole, 
black hole prevents the gas density from reaching the In the central 30 pc, two young dense star clusters, the 
Roche value high enough for star formation through self Arches and Quintuplet clusters, are observed (Nagata et 
gravity. To solve this problem, two possibilities have been al. 1995; Figer 2004). This model can naturally explain 
suggested: (1) star formation in situ in a dense gaseous the existence of two stellar disks. Numerical simulations 
disk around Sgr A*, or (2) infall of a young star cluster. of this model have been performed (Portegies Zwart et 
The disk model was proposed by Levin & Beloborodov al. 2003, hereafter PZ03; Kim & Morris 2003; Giirkan & 
(2003). A dense gaseous disk is formed, when a infalling Rasio 2005). These simulations showed that it took too 
molecular cloud is tidally disrupted by the central BH. If long for the star cluster to sink to the central parsec un- 
the density exceeds M/2'kR^ at radius R, where M is the less it was very massive (> IQ^Mq) or the initial position 
mass of the central BH, the disk becomes unstable with of the star cluster is very near from the GC (2pc or less). 
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In these simulations (PZ03; Giirkan & Rasio 2005), the 
dynamical friction fr-om the Galaxy on the star cluster 
was calculated analytically using the dynamical friction 
formula (Chandrasckhar 1943). However, it is not clear 
whether this analytical treatment of the orbital evolution 
is correct in the case of star clusters. For rigid objects, the 
orbital evolutions obtained using the analytical treatment 
agree well with those obtained using iV-body simulations, 
if In A is determined correctly taking into account the dis- 
tance from the galactic center and the size of the object 
(Hashimoto ct al. 2003; Spinnato et al. 2003). However, 
our fully self-consistent A^-body simulation of a satellite 
galaxy within its parent galaxy showed that the orbital 
decay of the satellite is much faster than those calculated 
analytically using the dynamical friction formula, even if 
the correct value of In A is used (Fujii et al. 2006). This 
difference is caused by the enhancement of the dynami- 
cal friction by particles escaped from the satellite. In the 
case of star clusters, the same enhancement should occur. 
Thus, a fully self-consistent A^-body simulation is neces- 
sary to obtain correct results for the orbital evolution of 
star clusters. 

Such a fully self-consistent A^-body simulation has not 
been performed because it was impossible with existing 
numerical methods. There are two classes of numerical 
methods for A-body simulations. One is the direct sum- 
mation, suitable for coUisional systems, and the other is 
schemes for collisionless systems such as tree, SCF, and 
particle-mesh schemes. Neither of them can handle such 
multi-scale systems like star clusters within galaxies. 

The direct summation is the simplest and the most ac- 
curate way to calculate forces; for a particle we calculate 
A^ — 1 inter-particlc forces and sum them up. Thus the 
cost to integrate all A^ particles is 0{N'^). The combi- 
nation of direct summation, Hermite integrator and indi- 
vidual timesteps is widely used. It can achieve very high 
accuracy. Hence, this method is used for simulations of 
collisional systems such as star clusters, planetary forma- 
tion, and galactic nuclei with central BHs. However, we 
cannot use this method for galaxies because its cost is too 
high. 

The tree algorithm (Barnes & Hut 1986; Makino 2004) 
is an approximate algorithm to calculate forces. It can 
achieve 0(A^log7V) scaling. This method is often used 
with 2nd-order leapfrog integrator. This method is very 
powerful for systems which need many particles and in 
which particles have similar timcscales, such as galaxies 
or large scale structures. Tree code is difficult to com- 
bine with individual timesteps and high order integration 
schemes. Therefore, the tree algorithm is not suitable for 
collisional systems. Other schemes for collisionless sys- 
tems have similar proper times. 

Fully self-consistent A-body simulations of star clusters 
which orbit in their parent galaxy require (1) accurate 
time integration for star clusters and (2) fast time inte- 
gration for the parent galaxy. Neither of direct or tree 
(or collisionless) methods, however, can satisfy these re- 
quirements. This is the reason why fully self-consistent 
A-body simulations of star clusters within galaxies have 



never been performed. 

A way to solve this problem is an algorithm based on 
tree code with individual timesteps (Hernquist & Katz 
1989; McMillan & Aarseth 1993). This approach, how- 
ever, has problems with the accuracies and/or the costs. 
Furthermore, it is difficult to use GRAPE with such a 
method. 

We developed a new algorithm which is an extension of 
the mixed variable symplectic (MVS) method developed 
for long-term integration of planetary systems (Wisdom & 
Holman 1991; Kinoshita et al. 1991). In our new scheme, 
star clusters are integrated using direct summation and 
Hermite scheme with individual timesteps, while galaxies 
are integrated using the tree code and leapfrog scheme 
with shared timestep. We combine them by embedding 
direct scheme into the tree code. Using this new algo- 
rithm, we have made it possible to perform fully self- 
consistent A'-body simulations of star clusters which orbit 
in their parent galaxies with very accurate time integra- 
tion for star clusters and fast time integration for the par- 
ent galaxy. 

We describe our new algorithm in section 3 after a de- 
scription of the MVS method in section 2. In section 4, 
we show the results of test simulations and performance 
of our new algorithm. We summarize in section 5. 

2. The Mixed Variable Symplectic Method 

The MVS integrator were introduced by Wisdom & 
Holman (1991) and by Kinoshita, Yoshida, & Nakai 
(1991). It is now widely used for long-term integrations of 
planetary systems. In the case of planetary systems, their 
Hamiltonian can be divided into Kepler motions and inter- 
action between planets. The MVS integrator suppresses 
the error of the motion due to the numerical integration of 
the solar potential since it integrates the Kepler motions 
analytically. 

Let us first describe briefly a symplectic integrator, the 
leapfrog scheme. The Hamilton equation is rewritten in 
terms of a Poisson bracket operator as 



(1) 



where / is a function of t. If we introduce a differential 
operator Dh defined as Dnf-— {f,H}, the formal solution 
of equation (1) is written as 



fit) 



(2) 



An integration algorithm can be thought of as an approx- 
imate expression of this operator. 

As an example, we describe a second-order leapfrog inte- 
grator. The Hamiltonian for an A^-body system is written 
as 

EPi \ ^ Gminij 
2m" 



(3) 



If we define 



No. ] 

H^ = -f2^^ (4) 

N 2 
i 

then we can express the the formal solution, the time evo- 
lution from < to i + Ai, as 

where A:— Dha ^-nd B:— Dhb- This operator can be 
rewritten by the Taylor series as 

fc 

^At(A+B) ^ -Qga.AtAgb.AtB ^ 0(Ar+l), (7) 

i=l 

where (ai^bi) (i — 1,2^ . . . ,k) is a set of real numbers and 
n is a given integer, which corresponds to the order of 
integrator. By neglecting the 0(A<"+^) then, we obtain 
a mapping from f{t) to f'{t + At) as 

fc 

/'(< + Ai) = ne'^'^*V'^*^/(0. (8) 

4=1 

This mapping is symplectic because it is just a product 
of symplectic mappings. This is an n-th order symplectic 
integrator. We can achieve n = 2 with k=2, with the choice 
of the coefficients ai — a2 — 1/2 and bi — l,b2 — 0- Now, 
equation (7) is reduced to 

^At(A+B) ^ giAtAgAtBgiAtA ^ 0{At^). (9) 

Therefore the time evolution is expressed as 

fit + At) = e^^*^e^*^e*^*^/(t). (10) 

This is the second-order leapfrog scheme, which is rewrit- 
ten as 

vi=vo + l- At ao, (11) 
^ z 

xi=xo + At vi, (12) 

vi=vo + ^ At ai, (13) 

where subscripts, 0, ^, 1, correspond to values at t,t + 
i At, t + At, respectively. 

The procedure of leapfrog scheme is as follows. 

1. Calculate the acceleration at a time, <, and update 
the velocity [eq. (11)]. 

2. Update positions using new velocity vi [eq. (12)]. 

3. Calculate the acceleration at t + At using the new 
positions, Xi, and update velocity [eq. (13)]. 

4. Repeat 1-3. 

Now, we explain an MVS integrator. The Hamiltonian 
for a planetary system can be expressed as 

H = HKcp + Hlnt, (14) 

where Hkcp is the kinetic energy plus solar potential and 
-ffint is the interaction energy between planets. If we define 
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Ha = Hint, Hb ~ Hkcp, (15) 
equation (10) becomes 

fit + At) = e^^"e^*^e^^*V(<)- (16) 

Here I-.—Dni^t and K:— Dh^^^p- Note that e^^^ generates 
motions of planets along unperturbed Kepler orbits, while 
gAt/ generates changes of momenta due to planet-planet 
interactions. This changes of momenta are called "velocity 
kicks." 

The difference from the usual leapfrog integrator is that 
pAtK gjygj]^ analytically by Kepler motion. Therefore 
the MVS method is expressed as 

v'l = vq + 1- At aint.o, (17) 

2 2 

a^o ~> (Kepler motion) Xi, (18) 
vi (Kepler motion) v'l, (19) 

vi = v\+l- At aint,i. (20) 

The integration proceeds as follows: 

1. Calculate the accelerations of planets due to gravi- 
tational interactions between planets, aint,0: at time 
t and change velocities by giving the velocity kicks 
[eq.(17)]. 

2. Update the positions and the velocities by At along 
its osculating Kepler orbit analytically [eq.(18) and 
(19)]. 

3. Calculate aint,o at t + At and change velocities by 
giving the velocity kicks [eq.(20)]. 

4. Repeat 1-3. 

MVS is a very powerful algorithm for long-term inte- 
gration of planetary systems. In general, the integration 
errors are 0{At"), where n is the order of the integrator. 
With a MVS integrator, if Hint is 0{e) of i^Kcpicr, the in- 
tegration errors are only 0(eAt"). This e is of the order of 
the planetary mass in the unit of solar mass and usually 
is very small. As a result, the error become much smaller 
than that of usual symplectic methods. 

3. The New Hybrid Scheme 

Now we consider simulations of systems consisting of a 
star cluster and its parent galaxy. For such a simulation, 
our new scheme should provide: 

1. high accuracy for star clusters 

2. fast integrator for galaxies 

3. fully self-consistent treatment of the total system. 

To achieve these goals, we constructed a new scheme 
which is a combination of the direct and the tree scheme. 
In our new scheme, the internal interactions of star clus- 
ters are calculated with high accuracy using the direct 
and Hermite scheme, while all other interactions (galaxy- 
galaxy, galaxy-star cluster) are calculated with the tree 
algorithm. We combine these two methods by extending 
the idea of the MVS. 
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In the MVS scheme, the Hamiltonian is divided into the 
kinetic energy plus solar potential and the interaction en- 
ergy between planets. In our hybrid scheme, we separate 
the Hamiltonian as 



H ~ Ha + Hp, 



Hb 



Ng 

y — 

2 — 1 ' 



Nsc 

E 

i=l 



PSG 



Ng Nsc ^ 

Ns 



(21) 
(22) 



2rnsc,; 



where iVc and A'^sc are the number of the galaxy particles 
star cluster particles, respectively, and Ha is the potential 
energy of the gravitational interactions between galaxy 
particles and between galaxy particles and star cluster 
particles, while Hp is the kinetic energy of all particles 
and the potential energy of star cluster particles. Using 



the replacement Ha 
time evolution as 



Ha and Hb = Hr, wc obtain the 



fit + At) 



giAtagAt/JgiAta 



fit), 



(24) 



where a:=DH^, P-= Dh^- 

In our scheme, we integrate star cluster particles and 
galaxy particles in different ways. Let us first discuss star 
cluster particles. They are integrated in a way similar to 
the MVS. The Keplerian in the MVS corresponds to the 
second and third terms in equation (23). Unlike the MVS, 
however, we cannot solve the Hamiltonian analytically. 
Hence, we replace an analytical solution (Kepler motion) 
in the MVS by a solution calculated by a higher-order 
integrator (e.g. the fourth-order Hcrmite integrator with 
individual timesteps). 

Thus, the integrator for star clusters is written as 



, 1 . 

xc,c,Q (Hermite scheme) x^ci, 



-'SCO 



(Hcrmite scheme) v 



SC,1' 



1 . 

- At a{G^sc4}' 



(25) 

(26) 
(27) 

(28) 



where subscripts, SC and G, stand for the star cluster 
and the galaxy, subscripts, 0, i, and 1, indicate times to, 



ti 



tf) + hAt, and ii = to + ^t, respectively, and v'^ 



sc,4 



represent a new velocity at ti, which have been integrated 
using the Hermite scheme. 

For galaxies, we use the leapfrog integrator expressed 

as 

1 
2 

1 . 

^G i + o a{All^G,l}, 



^G,i " '"G,0 + 77 At a{All^G,0}, 



Vg,1 



(29) 
(30) 

(31) 



where a^Aii^G} denotes the acceleration due to gravita- 
tional forces from all particles (including star cluster par- 
ticles) to the galaxy particle. The galaxy particles have 



longer timcscale than that of particles in the star clus- 
ter. Therefore, we adopt a second-order leapfrog integra- 
tor with shared timestep and tree algorithm. This scheme 
is less accurate than fourth-order Hermite scheme, but 
much faster and is symplectic. 

We call our new scheme "the Bridge scheme" (Bridge is 
for Realistic Interactions in Dense Galactic Environment). 
The procedure of the Bridge scheme is summarized in fig- 
ure 1 and as follows: 

1. Make a tree at to and calculate the accelerations 
from all particles on galaxy particles, a{Aii^G.o}i 
and from galaxy particles on star cluster particles, 
a{G^sc,o}, using the tree. 

2. Star cluster: give a velocity kick [eq. (25)]. 
Galaxy: update velocity [eq. (29)]. 

3. Star cluster: integrate positions and the velocities 
from to to ti using Hermite scheme with individual 
timesteps [eq. (26) and (27)]. 

Galaxy: Update position with leapfrog scheme [cq. 
(30)]. 

4. Make a new tree at ti and calculate the accelerations 
from all particles on galaxy particles, a.{Aii-^G.i}: 
and from galaxy particles on star cluster particles, 

a{G^SC4}- 

5. Star cluster: give a velocity kick [eq. (28)]. 
Galaxy: update velocity [eq. (31)]. 

As shown in equations (25) and (29), the forces on par- 
ticles of galaxies and those on particles of star clusters 
are calculated differently. The former is from all parti- 
cles, while the latter is from particles in the galaxy only. 
Therefore, we assigned two values of mass to each tree 
node. (We used center-of-mass approximation to forces 
from tree nodes.) One is total mass of all particles under 
the node, and the other is the mass of galaxy particles. 
To calculate forces on galaxy particles, we use the mass of 
all particles, and for star cluster particles we use the mass 
of galaxy particles. 



Leapfrog 



Free motion Free motion 
-.fc 



Velocity change 



Velocity change 



Velocity change 



MVS 



Kepler motion Kepler motion 



Perturbation 



Perturbation 



Perturbation 



Bridge 



Galaxy : Free motion 
Star cluster : Direct integration 



Free motion 
Direct integration 



Tree force 



Tree force 



Tree force 



Fig. 1. Procedure of the Bridge scheme. 
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4. Accuracy and Performance 

4.1- Comparison with Direct Scheme for Small-N Model 

As a test of the Bridge scheme, we performed a fully 
self-consistent A-body simulation of a star cluster within 
a galaxy and compared the results with those obtained 
with the fourth-order Hermite scheme. In this section, 
we describe the results and the performance of the Bridge 
scheme. 

We adopted a King model with the non-dimensional 
central potential Wq = 9 as the model of the parent galaxy 
and with VFo = 7 as that of the star cluster. The system of 
units is the Hcggic unit (Heggic & Mathicu 1986), where 
the gravitational constant G is 1 and the mass and the 
binding energy of the parent galaxy are 1 and 0.25, re- 
spectively. The initial position of the star cluster is at 
distance 2.5 from the center of the parent galaxy and the 
initial velocity is 0.65. Both the galaxy and the star clus- 
ter are expressed as a self-consistent A-body models. The 
number of particles of the parent galaxy, Ac is 10^ and 
that of the star cluster, Asc, is 2 x lO'^. In table 1, we 
summarize the model parameters and initial conditions. If 
we assume the total mass of the star cluster Msc = IQf'MQ 
and the unit length is 10 pc, the total mass of the Galaxy 
Mq = IO^Mq and the unit time and velocity are 0.15 Myr 
and 66 km s~^, respectively. These values would corre- 
spond to the central region of a galaxy somewhat smaller 
than our galaxy. 

The potential is softened using the usual Plummcr soft- 
ening. The softening length for the gravitational interac- 
tions between star cluster particles is esc = 2.0 x lO^'' — 
0.1 X 4/Asc, and that between others (i.e. between galaxy 
particles and for galaxy particles and star cluster parti- 
cles) is ec = 6.25 x lO'^. 

We used the opening angle = 0.75 with the center- 
of-mass (dipole-accurate) approximation. The maximum 
group size for GRAPE calculation (Makino 1991) is 8192. 
For the leapfrog integrator, we adopted the stepsizes of 
At = 1/128 and 1/256. The maximum timestep for the 
Hermite scheme with individual timesteps is equal to the 
timestep of the tree. All particles synchronize at each 
timestep for tree. Within these steps, star cluster par- 
ticles are integrated the Hermite scheme with individual 
timesteps (Makino & Aarseth 1992). For the timestep cri- 
terion, we adopted the standard formula, which is given 
in Makino & Aarseth (1992). These parameters are sum- 
marized in table 2. 

The simulations are performed on GRAPE-6 (Makino 
et al. 2003) for runs Direct 1 and 2, GRAPE-6 A 
(Fukushige et al. 2005) for runs Bridge 1 and 2. We sum- 
marize them in table 3. Total energy was conserved to 
better than 0.06% when M = 1/256, 0.2% At = 1/128 
with the Bridge scheme, and 0.006% with the Hermite 
scheme. 

To compare the results, we followed the time evolution 
of the position, bound mass, core radius, and core density 
of the star cluster. The bound mass and orbit are calcu- 
lated in the same way as in Fujii et al. (2006). We defined 
the core radius and the core density using the formula 



proposed by Casertano & Hut (1985). 



Table 1. Model Parameters of Testmodel 



Parameters 


Galaxy 


Star cluster 


Galactic halo 


King 9 


King 7 


Total mass 


1.0 


0.01 


Binding energy 


0.25 


2.5 X 10-" 


Half-mass radius 


9.8 X 10-1 


8.1 X 10-2 


N 


10^ 


2 X 10^ 


Initial position 




( 2.5, 0, 0) 


Initial velocity 




( 0, 0.65, 0) 



Table 2. Parameters for A'^-body Simulation 



Parameters 


Value 


ec 


6.25 X 10"^ 


esc 


2.0 X IQ-^ 


e 


0.75 


^crit 


8192 


At 


1/128, 1/256 



Figures 2 and 3 show the evolution of the distance from 
the galactic center and the bound mass of the star cluster. 
Figures 4 and 5 show the core radius and the core density 
of the star cluster. These results show that the Bridge 
scheme works very well. The difference between the re- 
sults of the Hermite scheme and that of the Bridge scheme 
is smaller than run-to-run variations in each method. 

Figures 4 and 5 show that core collapse occurs a,t T — 
150 — 180. Core collapse occurs at 



tr. 



'- Ctrh, 



(32) 



where trh is the star cluster's half-mass relaxation time 
(Spitzer & Hart 1971), 

Here rh, iVgc, and Afgc are the half- mass radius, the 
number of particles, and the mass of the star cluster. We 
adopted the Coulomb logarithm InA ~ In(O.lAgc). In 
isolated star cluster in which all stars has the same mass, 
c ~ 15 (Cohn 1980). From these equations, the core col- 
lapse time of our cluster is calculated as tec — 180. This 
value is consistent with the results of our simulation. Note 
that in our model we used the scaling of Mq = 4i?G = 1 ■ 

The total energy error of the system in Bridge 2 run 
is shown in Figure 6. The total energy is conserved very 

Table 3. Runs 



Runs 


methods 


seed 


stepsize 


run time (h) 


Direct 1 


direct 


1 


1/256 


34 


Direct 2 


direct 


2 


1/256 


34 


Bridge 1 


hybrid 


1 


1/128 


10 


Bridge 2 


hybrid 


2 


1/256 


19 
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Fig. 2. Distance of the star cluster from the GC plotted as 
a function of time. Solid and dashed curves show the re- 
sults of the runs in which all particles calculated with the 
direct (Hermite) scheme. Dashed-doted and dotted curves 
show those with the Bridge scheme. 
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Fig. 3. Bound mass of the star cluster plotted as a function 
of time. Curves have the same meanings as in figure 2. 



Fig. 4. Core radius, rc, plotted as a function of time. Curves 
have the same meanings as in figure 2. 
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Fig. 5. Core density, pc, as a function of time. Curves have 
the same meanings as in figure 2. 

well. The total energy error depends on the parameters 
for tree, Ai and 0. This is because most of the energy 
error is generated in the parent galaxy, which is much 
larger than the star cluster and has much larger energy. 

To see whether the internal energy of the star cluster 
is conserved or not, we measured the energy error of the 
internal motion of the star cluster within each step, Af. 
The cumulative error of each step is shown in Figure 6. 
Note that wc plot the energy error of the star cluster rela- 
tive to the internal energy of the star cluster, which is 0.1 
% of the total energy of the system. Although the error 
become larger after core collapse occurred, it is conserved 
well. 

The distributions of the CPU time in runs with At = 
1/128 and At = 1/256 are shown in table 4. In these 
simulations, the cost of the tree part was much larger than 
that of the direct part. The CPU time of the tree part is 
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'o r 1 

-" 50 100 150 200 

T 

Fig. 6. The energy error of the system, for the hybrid run 
with At = 1/256 (Run Bridge 2). Full curve shows the to- 
tal energy error of the system and dashed curve shows the 
cumulative energy error of the star cluster. 

almost constant throughout the simulation. In contrast, 
the cost of the direct part increased after T ~ 140, because 
the core density become higher. 



Table 4. Distribution of the CPU Time of our Simulations 


Section of Code 


Percentage 


of CPU Time (%) 




At 1/128 


At = 1/256 


Tree 


87.2 


91.2 


Direct 


10.7 


6.6 


The others 


2.1 


2.2 



4.2. Large-N Models 

We performed simulations of large- TV models. The num- 
ber of particles is sufficiently large for simulations of star 
clusters near the GC (iVc = 2 x 10^, N^c = 65,536). The 
model of the galaxy represents the central region of the 
Galaxy from ~ 1 to several pc from the GC. For the star 
cluster, we followed the model used in Portegies Zwart et 
al. (2003). They modeled the Arches and Quintuplet star 
clusters. The core radius of our model, rcorci is 0.087 pc. 
Using this model, we performed two simulations, in which 
the star cluster has circular and eccentric orbits, respec- 
tively. In table 4.2, we summarize the model parameters. 

We performed iV-body simulations using the Bridge 
code. For the tree part, we used the opening angle 9 = 
0.75 with the center-of-mass (dipole-accurate) approxima- 
tion. The maximum group size for a GRAPE calculation 
(Makino 1991) is 8192. The stepsize of leapfrog integra- 
tor is At = 1/512 (Heggie unit). The potential is softened 
using Plummer softening. The softening length for gravi- 
tational interactions between star cluster particles, egc, is 
1.0 X 10^^ pc and that for others, ec, is 3.9 x 10^^ pc. We 
stopped the simulations at T — 0.75(Myr) = 5(unit time). 
These parameters are summarized in table 6. After the 



core collapse, the structures of the star clusters are not 
expressed correctly in our simulations because we use a 
softened potential for stars. 

We used GRAPE-6 (Makino et al. 2003) for force cal- 
culation. The total energy was conserved better than 
5 x 10~^ for the circular orbit (figure 9) and 8 x 10~^ for 
the eccentric orbit (figure 10) throughout the simulations. 



Table 6. Parameters for A'^-body Simulation 



Parameters 


Value 


ec 


3.9 X 10-^ (pc) 


esc 


1.0 X 10-5 (pc) 


At 


2.9 X 10-"* (Myr) 


(Heggie unit) 


1/512 


6 


0.75 


«crit 


8192 



Table 7. Initial Conditions 



Simulation Initial position (pc) Initial velocity (km s 
Circular 2 130 

Eccentric 5 72 



Figure 7 show the snapshots from the run in which the 
orbit of the star cluster is eccentric. Figure 8 shows the 
time evolution of the distance from the GC, bound mass, 
and core radius of the star clusters. In both simulations, 
core collapse occurs at 0.5 - 0.6 Myr. We obtained the 
core collapse time, tec = 0.51 Myr, from equation (32) 
and (33), where the half-mass radius of the star cluster, 
rh = 0.13 (pc). We adopted c = 0.20, which is suggested 
by (Portegies Zwart & McMillan 2002). The core coUapse 
time of our simulations is consistent with the results of 
the previous studies. 

Figure 9 and 10 show the total energy error of the sys- 
tem and the internal energy error of the star clusters. The 
internal energy errors are cumulative as in small- mod- 
els. The energies conserve very well. 

The total CPU times and the distributions of the CPU 
time are shown in table 8. The total CPU time was about 
40 hours. The direct part consumes about half of the CPU 
time. 

4-3. Performance Model of the Hybrid Scheme 

We analyzed the CPU time of each part in detail. 
Figure 11 shows the CPU time per 4 steps for each part for 
the run with circular orbit. Simulation time is represented 
using the Heggie unit. We used the parent galaxy to de- 
fine the Heggie unit, i.e., Mq = 4Eq = 1. The unit time 
in the Heggie unit corresponds to 0.15 Myr. Hereafter we 
use the Heggie unit for time to discuss the performance of 
the hybrid scheme. The CPU time of the tree part is al- 
most constant throughout the simulation. In contrast, the 
cost of the direct part gradually decreases and suddenly 
increase after T ~ 3.5. 
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Fig. 7. Snapshots of the star clusters projected onto x — y plane. The orbit of the star cluster is eccentric. 
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0.66 
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Fig. 8. The distance from the GC (top), bound mass (middle), and core radius (bottom) of the star clusters plotted as a function 
of time. The orbits of the star clusters are circular in the left panels and eccentric in the right panels. 



Table 8. CPU percentage for the test models 



Section of Code 


CPU Time (sec) 
Circular Eccentric 


Percentage of CPU Time (%) 
Circular Eccentric 


Tree 
Direct 
The others 


5.6x10'' 6.3 xlO'' 
7.5 X 10" 7.5 X 10^ 
1.1 xlO^ 1.4 xlO^ 


42.3 45.3 
56.9 53.6 
0.8 1.0 


Total 


1.3 X 10^ ~ 37(h) 1.4 x 10^ -39(h) 


100.0 99.9 



10 



Author(s) in page-head 




o 



o 
o 



o 

in 



o 
o 



o 



1 




: Circular orbit 


Direct . 


- 

- 


I'l T 

|! «fj' * 

f : 




Tree- 








The others - 



T 



T 



Fig. 9. Same as figure 6, but for large- A' model with circular 
orbit. 
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Fig. 10. Same as figure 6, but for large- Af model with eccen- 
tric orbit. 



As shown in figure 13, the CPU time of the direct part 
is proportional to the number of the steps for the Hermite 
scheme, rtstop- Figure 12 shows the time evolution of the 
average number of timesteps per particle, rigtep- It gradu- 
ally decrease until T = 3, and stays nearly constant. In fig- 
ure 13, CPU time suddenly increases starting at T = 3.6. 
This time corresponds to the time of core collapse, and 
after that the internal dynamics of the star cluster is not 
correctly followed because of the finite softening. In the 
discussions below, we consider the behavior of the CPU 
time before core collapse. 

From these results, we can construct the performance 
model of the Bridge code. The total CPU time for a step, 
Ai, can be written by the number of the tree particles, 
A^tieej that of the direct particles, A'diroctj a-nd that of steps 
for Hermite scheme per step, ristcp- The cost of tree is 



Fig. 11. CPU time of the direct and tree parts per 4 steps 
(At = 1/128) for the circular orbit. The solid, dashed, and 
dotted curves show the CPU time of tree, direct, and the 
others parts, respectively. 




T 



Fig. 12. The number of steps of the direct part per particle 
per 4 steps (1/128 unit time) for the circular orbit. 
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Fig. 13. The CPU time of the direct part per 4 steps (1/128 
unit time) before T = 3.6. Filled circle and open circle show 
the results of the circular orbit and eccentric orbit. Solid line 
shows the model of equation (34). 

proportional to A'tree log iVtrco ^ Mice- The CPU time of 
the direct part depends on A'trco and nstep. The cost of 
force calculation is proportional to A^liroct ^^'i the other 
costs are proportional to iVdircct- Therefore, the total CPU 
time per At is given by 



■step 7 

(34) 

where a, /3, and 7 are constants. Here, a is almost con- 
stant through a simulation, but depends on 9. The value 
of 7 is determined by the performance of GRAPE. Makino 
et al. (2003) shows that the calculation time on GRAPE 
per interaction per particle is expressed as 

Tgrape = » , , , „7__ (sec), (35) 



9 X lO^r 



'pipes 



where 



''pipes 



is the total number of pipelines. With 



GRAPE-6A ripipcs = 24, with GRAPE6 ripipes = 192. 
Hence, we estimated 7 = 4.6 x 10~^° (sec), for GRAPE- 
6A, 7 = 5.8 X 10"" (sec) for GRAPE-6. From the re- 
sults of our runs, we estimated the values of the constants, 
a = 1.2 X 10-5 (sec) for 9 = 0.75 and /3 = 6.9 x 10"^ (sec). 

The number of particles that we need for fully self- 
consistent A^-body star cluster simulation is Nq ~ 2 x 10^ 
for a galaxy and N^c ^ 6.5 x 10^ for a star cluster. In this 
case, the total CPU time for the Bridge scheme is esti- 
mated as 1.2 X 10^ sec ^ 34 hours for a simulation with 
At = 1/512 and 5 unit time integration on GRAPE-6. 
Here we used Ugtcp ~ 33 from the results of our simula- 
tions. The actual time for such a simulation was 37 - 39 
hours (see table 8). So our model predicts the CPU time 
with ~ 20 % accuracy. 

We can also estimate that for the Hermite scheme. The 
CPU time of the Hermite scheme, or the direct scheme, 
can be estimated using the second and third terms of 
equation (34), where we used ngtep = 770 per unit time. 
Therefore, the total CPU time is estimated as about 260 



hours per 5 unit times on GRAPE-6 for iV = 2 x 10''. It's 
about seven times longer than that for the Bridge scheme. 

5. Summary and Discussion 

5.1. Summary 

We have developed a fast and accurate algorithm, "the 
Bridge scheme," for fully self-consistent A^-body simula- 
tions of a star cluster moving in its parent galaxy, where 
both are modeled as A^-body systems. The Bridge scheme 
is a hybrid of the tree and direct schemes and is based on 
an extension of MVS. Wc performed self-consistent A^- 
body simulations of a star cluster in a galaxy and com- 
pared the results with the Bridge scheme and that with 
the direct scheme (the Hermite scheme). They agreed 
each other very well and the energy error was sufficiently 
small. We also showed that we can perform a full A^- 
body simulation of a star cluster and a galaxy system with 
A^SC = 65536 and A"g = 2 x 10^ using our new scheme more 
than seven times faster than the direct scheme. 

5.2. Comparison with Tree-based Algorithms 

In previous studies, several tree-based algorithm with 
block timesteps were developed. Hernquist & Katz (1989) 
adopted block timesteps in a tree code. In their scheme, 
the tree is reconstructed in each step. When the timesteps 
of particle do not vary so widely, the cost of the tree re- 
construction is not so expensive. However, the cost is very 
expensive for star clusters, because star clusters have wide 
range in their timesteps. 

McMillan & Aarseth (1993) developed a tree-based 
high-order integration scheme for collisional systems using 
block timesteps and multipole (up to octupole) expansion. 
In this scheme, the tree is reconstructed at the appropriate 
cell timesteps determined by the motions of the particles 
in the cells. Instead of reconstructing tree at each step, 
the moment of each cell is predicted. However, the accu- 
racy is limited by the time interval of tree construction. If 
longer timesteps are permitted, the tree evolves during the 
step and the errors increase. If the time interval is short, 
the cost of tree construction become large. In addition, 
their algorithms are difficult to use with GRAPE. 

5.3. Applications to Other Problems 

Our initial motivation for developing the Bridge scheme 
is to use it for the problem of a star cluster orbiting in its 
parent galaxy. However, it might have much wide ap- 
plication range. For example, if the parent galaxy has 
the central massive black hole, it is natural to handle it 
and stars near by with the direct scheme, and the rest 
of the system by tree. In this case, some particles must 
move "tree" and "direct" treatment, but in principle such 
a code can be developed. Our method can be applied to 
any large- A^ systems in which small part of the system 
shows collisional behavior. 

The authors thanks Pict Hut for useful comments and 
the name of the hybrid scheme, Keigo Nitadori and 
Ataru Tanikawa for fruitful discussions, and the referee. 
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